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Abstract 


The  development  of  a  comprehensive  computational  fluid  dynamics  approach  for  conducting 
simulations  of  projectile  penetration  into  water  and  dry  sand  is  reported.  High  resolution  upwind 
schemes  suitable  for  a  fluid  dynamic  system  consisting  of  gas,  liquid,  and  dispersed  solids  phases  arc 
derived  and  are  combined  with  a  time-derivative  preconditioning  strategy  for  efficient  time  integration 
at  all  flow  speeds.  A  solids-stress  model  based  on  Mohr-Coulomb  critical-state  theory  is  used  to 
account  for  compaction  and  deformation  of  sand  during  projectile  penetration.  An  overset-mesh 
framework  is  also  implemented  in  order  to  handle  projectile  relative  motion  in  subsequent  work,  and 
improved  phase  interface  capturing  methods  are  also  developed  and  tested.  Results  arc  presented  for 
two  sets  of  experimental  data  involving  projectile  penetration  into  dry  sand.  The  computational  results 
are  sensitive  to  the  solids-stress  model  and  the  drag  coefficient  predictions  are  generally  lower  than 
indicated  in  the  experimental  data. 
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1.  Background  and  Objectives 

A  problem  of  great  interest  for  naval  applications  is  the  water  and/or  sand  entry  problem  for  high¬ 
speed  projectiles.  Supercavitating  projectiles  can  be  used  for  underwater  mine  neutralization,  beach 
and  surf  zone  mine  clearance,  littoral  ASW,  and  neutralizing  combat  swimmers.  Many  current  and 
possible  future  systems  must  deal  with  a  high-speed  water  entry:  RAMICS,  VENOM,  HYDRA,  SEA 
SNIPER,  JABS,  and  anti-combat  swimmer  systems.  The  water  entry  phase  of  flight  is  interesting  and 
challenging  due  to  projectile  transitioning  from  flight  in  air  to  supercavitating  flight  through  water. 
For  the  beach/surf  zone  clearance  application,  projectiles/bombs  are  designed  to  enter  the  water  and 
sand  vertically,  so  that  an  axisymmetric  approximation  to  the  flowfield  may  be  appropriate  if  the 
projectile  is  stable.  Investigating  projectile  stability  will  require  a  fully  three-dimensional  approach 
with  a  six-degree  of  freedom  (6-DOF)  algorithm.  Projectile  stability  should  be  a  strong  function  of 
supcrcavity  size,  shape,  and  unsteadiness.  The  supercavity  is  largely  determined  by  velocity  and  nose 
shap.  The  nose  shape  of  a  penetrator  is  generally  a  compromise  between  lethality  and  cavity 
generation  concerns,  with  minimizing  drag  in  air  being  a  tertiary  consideration. 

The  overall  goal  of  the  presented  work  is  to  develop  the  ability  to  predict  velocity  and  penetration 
depths  in  water  and/or  sand  of  a  penetrator  or  bomb  after  a  vertical  water  entry.  The  work  seeks  to 
develop  new  algorithms  for  simulating  the  water  and  sand  entry  phases  of  high-speed  projectiles.  The 
ability  to  accurately  predict  velocity  and  penetration  depth  is  critical  to  lethality  of  systems  that  require 
a  maximum/minimum  velocity  on  impact  or  a  specific  sand  penetration  depth.  The  current  work 
develops  high-resolution  upwinding  schemes  and  time-derivative  preconditioning  techniques  suitable 
for  a  fluid-dynamic  system  consisting  of  gas  (air  or  cavitatcd  water),  liquid  (water),  and  dispersed  solid 
(sand)  phases.  Such  techniques  will  facilitate  sharp  capturing  of  phase  interfaces  and  vortical 
structures  at  all  flow  speeds  and  under  all  states  of  compressibility,  thus  enabling  high-fidelity 
predictions  of  the  forces  and  moments  acting  upon  high-speed  projectiles  under  transient  conditions. 

The  technical  approach  combines  the  air  /  water  two-phase  flow  Navier-Stokes  codes  developed  by 
NSWC-PC  and  NCSU  [1-6]  with  the  dense  gas  /  solid  fluidization  codes  developed  at  NCSU  [7,8]  to 
yield  a  unified  procedure  capable  of  capturing  the  detailed  physics  of  projectile  penetration  into  water 
and  into  dry  sand.  The  complete  algorithm  allows  detailed  modeling  of  granular  flow  effects,  such  as 
solids  stresses,  and  captures  phase  interfaces  as  part  of  the  solution.  Velocity  slip  effects,  solids 
compaction,  and  cavity  bubble  formation  and  collapse  are  accounted  for  in  the  model.  The  new 
algorithm  utilizes  time-derivative  preconditioning  as  a  way  of  reconciling  the  widely-varying 
characteristic  speeds  found  in  this  three-phase  system,  thus  enabling  efficient  integration  at  all  flow 
speeds  and  all  states  of  compressibility.  To  enable  eventual  application  of  the  code  in  six-degree  of 
freedom  simulations,  the  procedure  is  embedded  within  an  overset-mesh  framework  using  the 
SUGGAR/DiRTLIB  [9,10]  protocol.  The  sections  that  follow  describe  the  governing  equations,  the 
numerical  methods  used  to  solve  them,  and  the  results  of  this  investigation,  followed  by  some 
concluding  remarks. 
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2.  Governing  Equations 


The  governing  equations  for  the  three-phase  system  (fluid  (air  +  water)  and  sand)  may  be  written  in 
tensor  notation  as  follows 


Fluid  phase: 

vapor  mass: 


fluid  mass: 


C  Y  ?P  !  d(PfafK)  t  diPfafY^jj)  Q 
J  v  dr  dt  dxj 

c  dp  |  8(pfaf)  ,  djPfCtfUjj )  Q 

cr  dt  8xj 


fluid  momentum:  0fui  ,  —  + 

(  ,J  5t 


dp  djpj^u^/)  dJ^ctjU^jUjj-T^ 


dt 


dx , 


=  -a, 


dp_ 

dx. 


+  prargi+S„ 


fluid  energy: 


0fH  f  —  + 

f  f  8r 


dp  d(pfafHf)  dp  t  d(pfafHf-T)k  fukf-qf)  _ 


dt 


dt 


dx. 


PfafUi.fgi  “Energy 


(1) 

(2) 

(3) 

(4) 


Solids  phase: 

solids  mass: 

solids  momentum: 


e  dPs  !  !  d(Psasuj,s)  Q 

s  dr  dt  dx  j 


Ps“s 


<Ks 

dt 


8U: 


dXjJ 


dp  dTy 

=  -as  —  +  —  +  Psasgi  -  Sn 
OX:  OX  , 


(5) 

(6) 


solids  energy:  psasCp 


'dT  dul  s  ^ 


=  s. 


energy 


(7) 


The  subscripts  v,f,  and  s  represent  vapor,  fluid,  and  solids,  respectively.  The  density  is  denoted  by  p , 
the  velocity  in  coordinate  direction  i  by  w, ,  the  volume  fraction  by  a ,  the  vapor  mass  fraction  (relative 
to  the  mass  of  the  fluid  phase)  by  Yv ,  the  pressure  by  p  ,  the  fluid  and  solid  temperatures  by  T  and  Ts , 
and  the  gravitational  acceleration  vector  by  g, .  The  terms  involving  %T  are  associated  with  the  time- 
derivative  preconditioning  method  discussed  in  Section  3.  The  system  is  solved  for 
V  =  [Yv,p,ui  f,T,as,ujs,Ts] .  Other  variables  arc  defined  in  the  next  section. 


2.1.  Fluid  Closure  Relations 


Closure  of  the  above  system  is  achieved  by  applying  suitable  equations  of  state  for  the  three  phases. 
For  the  fluid  phase,  the  following  thermodynamic  relations  hold: 


Amagat’s  law:  — 

Pf 


i-y„ 


Pv(pJ)  p,(pJ ) 


(8) 
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Mixture  total  enthalpy:  H f  =  Yvhv(T)  +  Ylhl(p,T)  +  ^(u 2  +v2  +  w2) 
Vapor  phase  equations  of  state:  pv  =  p/(RvT),  hv=hv(T)[  11] 


(9) 

(10) 


Liquid  water  equations  of  state  (modified  Tait  equation): 

p, = -  P,„(T))r, 


h,=h,(T)  + 


3  x10s 
(p  - 101325) 

Pi 


(ID 


Saturation-state  values  (subscript  ‘sat’)  are  obtained  from  [12]  and  the  liquid  reference  enthalpy  is 
obtained  from  [11], 


The  fluid  stress  tensor  is  modeled  as 


1  tokj  y 

2  dxkJ  ij 


Tijj=2MfafSiJJ , 

(  tov  ,  8uJ  f  ) 

(dXjj  dx,  f  j 

(12) 


with  the  fluid  viscosity  formulated  as 

Mf  =  avPv(T)  +  a,p,{T) .  (13) 

Here,  av  is  the  volume  fraction  of  the  vapor,  relative  to  the  fluid  volume.  The  vapor  and  liquid 
viscosities  arc  obtained  from  [11]  and  [13],  respectively. 

2.2.  Sand  Closure  Relations 


For  the  sand  (granular)  phase,  the  intrinsic  sand  density  ps is  taken  to  be  2676.605  kg/m  The  solids 
stress  tensor  77  is  defined  as 

Tf—pA+f&,’  <l4> 

where  the  solids  viscosity  jus  is  expressed  as  a  combination  of  a  Newtonian  viscosity  and  a  non- 

Newtonian  Mohr-Coulomb  [14-17]  viscosity  that  requires  that  the  frictional  stresses  exhibit  a  zero- 
order  dependence  on  the  fluid  strain  rate.  Considering  the  solids  pressure  as  the  critical-state  pressure 
leads  to  the  requirement  that 
Ps  sin(^) 


Ps=Ps.N  +- 


(15) 


JS:S\,  +£/</2  ’ 

where  (f>\s  the  angle  of  internal  friction  (taken  as  zero  or  28  degrees,  depending  on  whether  solids 


shear  stresses  arc  included.  The  strain-rate  tensor  for  the  solids  phase  is  given  as 


( 


So.s  -  2 


du. ,  du , 


\dXJ.s 


dx. 


du 


k,s 


i,s  y 


2&*,  " 


and  the  second  invariant  of  the  strain-rate  tensor  is 

S  ■  S  \s  =  Sl2,sS\2,s  +  +  ^2i,s^2J,s  ~  ^U,s^22.s  ~  S\  ~  ^22.s^3i.s 


(16) 

(17) 


6 


The  term  e Id2  (e  =  1  x  10  4,  d  p  =  0.00240  m)  is  used  to  prevent  divide-by-zero  in  uniform  regions  of 

the  flow.  The  solids-stress  model  can  be  viewed  as  the  distribution  of  the  solids  pressure  in  different 
directions  depending  on  the  elements  of  Sljs .  Thus,  it  is  appropriate  to  enforce  the  constraint  that 
S 

—  1  <— j= - !£ - -<1  (18) 

yjS'.Sl  +e/d2p 

for  the  non-Newtonian  part  of  the  stress  tensor.  As  Mohr-Coulomb  -  type  granular-flow  models  arc 
known  to  be  intrinsically  unstable  at  short  wavelengths  [18],  it  is  necessary  to  include  a  Newtonian 
viscosity  to  regularize  the  solution  and  to  prevent  unphysical  behavior.  We  have  adopted  a  simple 
model  based  on  a  von  Neumann  stability  analysis: 

Ms.n  =  c*.n  sinW Ps  <*sas  A ,  ( 1 9) 

where  A  is  a  characteristic  mesh-cell  dimension  in  a  particular  coordinate  direction,  Cs  N  is  a  model 


constant,  and  as  is  a  sound  speed  associated  with  the  granular  phase.  This  sound  speed  is  defined  as 


J_jK 

Ps  das 


(20) 


Results  presented  later  illustrate  the  effect  of  varying  the  model  constant  Cs  N  on  the  predictions. 


The  solids  pressure  /?vmust  be  defined  in  order  to  complete  the  formulation  of  the  solids  stress  tensor. 
The  solids  pressure  is  near  zero  when  the  particles  are  in  a  dilute  state  but  increases  rapidly  to  giga- 
Pascal  levels  as  the  sand  particles  are  compacted,  then  fragmented,  by  the  applied  load.  We  have  used 
two  techniques  to  model  the  response  of  dry  sand  considered  in  Lockheed-Martin’s  experiments  in  this 
study.  In  the  first,  we  use  an  approach  based  on  the  Coopcr-Eaton  equation  for  powder  compaction. 
[19]  The  Coopcr-Eaton  equation  relates  the  solids  pressure  (identified  as  the  trace  of  the  solids  stress 
tensor)  to  the  volume  fraction  and  has  the  general  form 

~  _ — ^o_  _  ^  exp(_^  / Ps )  +  a2  exp (-B2  / ps),  (21) 

^5. max  —  «,,0 

where  as  is  the  solids  volume  fraction  and  ps  is  the  solids  pressure.  The  Coopcr-Eaton  equation  is 


valid  for  solids  volume  fraction  values  greater  than  an  initial  unloaded  value  a,  0 .  As  written  above, 

the  Coopcr-Eaton  equation  is  implicit  in  the  solids  pressure.  To  avoid  having  to  solve  a  non-linear 
equation,  we  have  re-written  the  model  as 

a.  -  a , 


a,  -  -  A,  cm-B, rT^"] 

^.max-O'x.o  exp(-5, 1 ps) 

and  have  curve-fitted  the  term  in  brackets  as  a  fourth-order  polynomial: 
A2  exp(-5,  /  ps) 


(22) 


g(«5  )  =  [!+- 


-]  =  1,  as  <  0.4 


A i  exp(-5,  /  ps)' 

.  r,  A-,  exp(-#,  /  0.)n  _  ~2  ~3  ~4  ~  r.  X 

g(av)  =  [l  +  -2 -  2^]  =  a,  +a2as  +a}as  +aAas  +a,as,a5>0.4 

Ax  exp(-5,  / ps) 

This  enables  an  explicit  solution  for  the  solids  pressure: 


Ps.C-e,  I  —  " 


B , 


a. 


(23) 


(24) 


ln(— £— )-ln(4,) 
g(as) 
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To  extend  the  Cooper-Eaton  equation  smoothly  beyond  the  assumed  loading  condition  as  -  as  0 ,  the 

following  combined  form  is  used: 

„  ,  ,  .  dln(Ps.c-e,i  la,  =0.15)  ~ 

Ps.C-e. 2  =  eXP(ln(/>S.C-e.l  la,  =0.15) - - (°- 1  5  “  X  <*s  <  0. 1  5 


da. 


(25) 


Ps.C-e. 2  ~  Ps.C-e.\’as  ^0.15 

For  volume-fraction  values  less  than  a%  0 ,  we  the  Boivin,  et  al.  solids  pressure  model: 


Ps.b  =  PsCs,B[as  +  2a(  max  g  ln(l - %—)-  ] 

^s.max.B  ^s.max.B 


and  then  blend  between  the  two  to  yield  the  final  form  for  the  solids  pressure: 
Ps  =  exP{ln(ps.s)  +  -^[1  +  tanh(100(«5  -  av,0))](ln(pJ>c_£,2)  -  ln(pj  g))} 


(26) 

(27) 


Values  for  the  constants  in  the  Cooper-Eaton  model  were  initially  taken  from  data  given  for  silicon 
dioxide  in  [19]  (Model  I)  and  were  then  least-squares  re-fitted  to  match  sand  compaction  data  obtained 
by  Lockhecd-Martin  [20]  (Model  II).  These  values  are  shown  in  Table  1,  and  the  resulting  composite 
solids-pressure  profiles  are  shown  in  Figure  1 . 


Table  1:  Parameters  for  Cooper-Eaton  Solids-Pressure  Models 


Parameter 

Model  I  value 

Model  II  value 

as.O 

varies  depending  on  the 
case 

varies  depending  on  the 
case 

^5,  max 

0.93 

0.93 

A 

0.6 

0.6528 

A 

17.232  (MPa) 

67.033  (MPa) 

a\ 

-4.1175 

-4.1175 

a2 

32.8372 

32.8372 

ai 

-96.1131 

-96.1131 

aA 

121.4415 

121.4415 

a5 

-55.1670 

-55.1670 

CS.B 

1.397 

1.397 

&s,  max.  B 

1.0 

1.0 

In  some  of  the  later  calculations,  we  have  used  another  fitting  of  the  Lockhccd-Martin  sand 
compaction  data  in  an  attempt  to  provide  a  sharper  transition  between  loaded  /  unloaded  states  and 
thereby  improve  drag  predictions.  This  approach  is  defined  by  the  following: 


Ps  =CM2{AmK—)Bu'  -l]+^2[(— <*,><** 

a  si  a,, 

Ps  =  0.0,  as<as, 

and  the  constants  are  fitted  as  follows: 


(28) 
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Tabic  2:  Parameters  for  Refitted  Solids-Prcssure  Model 


Parameter 

Model  III  value 

°sl 

0.580769002438 

C M2 

6892700 (Pa) 

Am\ 

4.3634398407e-5 

Bm 

38.072977270364 

AM2 

4.41475836547218 

&M2 

7.14454835369866 

Solids  pressure  predictions  from  this  model  are  shown  in  Figure  2. 

2.3.  Interchange  Terms 

Momentum  transfer  between  the  fluid  and  solids  phases  is  modeled  using  classical  Stokes  /  Ergun  drag 
laws: 

Smom,  =  P,a,CA rag  («/,,  ~  t/„  )  ,  (29) 

For  af  >0.8,  a  Stokes  model  is  used: 

Cdrag=-C^.^-,  (30) 

y.s  “  p 

with 


Qn  .  - 


Pfafd> 


[1  +  0.15(2 


Pfafdp\v /~KtI\0.687  -i  P,<*f<*P\v,  -V,\ 


Pf 


bo/  j 


Pf 


<1000 


(31) 


Cdrv  =  0.44  \Vf-Vs  |,  ^^M>1000 
For  or,  <0.8,  an  Ergun  model  is  used: 

1 


^drag 


A«/ 


a  Mr  \Vf~V, 

150  — —2- +  1.75/?  .■  — - - 

t/J  7  dp 


(32) 


In  these  expressions,  the  relative  velocity  magnitude  is  given  as  |  Vf  -  Vs  |=  yj(ufJ  -m5  ,)(u,(  -  u~) . 


Energy  transfer  between  the  fluid  and  solids  phases  is  modeled  by  using  a  simple  relaxation  approach: 


PsaSCp.s 


du, 


8TS 

dt  +Ujy8x , 


=  S^rm  =  ~j2~kfNu  (T-Tp) 

dP 


energy  j2  J  « '  P> 

In  this  equation,  the  specific  heat  for  sand  (J/kg-K)  is  curve-fitted  as 

Cp.s  =as\  +aslTp  +asiTp  +asJ3p  +asSTp’ 

asl  = -1. 367634  xlO2 


(33) 


as2  =  3.950952 
as}  =  -4.94 1694x1  O'3 
a,4  =  2.868886x1  O’6 
as5  =  -6.1 8841 7x10-'° 


(34) 
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This  curve  fit  is  valid  for  temperatures  between  300  and  1600  K.  For  temperatures  below  300  K,  the 
specific  heat  is  fixed  as  674.7  J/kg-K,  and  for  temperatures  above  1600  K,  the  specific  heat  is  fixed  as 
1230  J/kg-K.  In  Eq.  (33),  kt  is  the  thermal  conductivity  of  the  fluid  phase,  and  the  Nussclt  number  is 

expressed  as  a  function  of  the  solids  voidage  and  relative  velocity  as 


Nu  =  (7  - 1 0a ,  +  5a) )(1  +  0.7  Re0  2  Pr1/3)  +  (1.33-  2.4a,  + 1 .2a2, )  Re0  7  Pr1/3 , 


Re  = 


Pf 


W-VA  d. 


Pr  =  0.72 


(35) 


Pf 


3.  Numerical  Methods 

3.1.  Eigenvalue  Analysis 

To  arrive  at  an  efficient  method  for  solving  the  system  of  equations  shown  above,  we  reformulate  the 
inviscid  components  of  the  equation  set  in  a  non-conservative  form  and  then  add  artificial  time 
derivatives  of  fluid  and  solids  pressure  to  the  continuity  equations  for  each  phase.  Writing  this  system 
in  a  quasi-linear,  one-dimensional  form  for  convenience,  we  have  the  following: 

(36) 

(37) 


5: 

9*1^ 

II 

O 

where 

VT  = 

[Yv,af,us,uf 

,P,T] 

and 

'  Pfaf 

0 

0  0 

( 

0 

-  Ps  (1  +  @sas  ) 

0  0 

( 

0 

0 

Psas  0 

( 

M  = 

0 

0 

0  pfaf 

( 

afpf\K 

pf 

0  0 

af(Pf\ 

Pfa/hf\y 

0 

0  pfafuf  pfafh) 

PfClfUf 

0  0  0 

0 

0 

~Ps“s  Psas  0 

0 

0 

~  Psa\  PsasUs  0 

0 

A  = 

0 

0  0  pfafuf 

af 

afU/Pf\y 

pfUf  0  Pfaf 

afufpf\h 

\ 

PfafUfhf  |, 

0  0  p  i'  a ,  m  y 

PfdfUfhf 

0 

0 

0 

0 

afpf\T 

pfafhfl 


0 

0 

0 

0 


afuf 


Pf\r 
Pfafufhfl 


(38) 


The  proper  specification  of  the  parameters  0f  and  ds  should  allow  control  over  the  variation  of  the 

eigenvalues  (characteristic  speeds)  of  the  system.  The  inclusion  of  the  solids  energy  equation  docs  not 
change  the  eigenvalues  and  thus  it  is  omitted  from  the  analysis.  Also,  it  should  be  noted  that  the 

r)n 

inclusion  of  the  fluid  pressure  -  a,  —  term  in  the  solids  momentum  equation  renders  the  system  non- 

ax. 
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hyperbolic  in  time  and  greatly  complicates  the  analysis.  As  such,  this  term  is  considered  to  be  a  source 
term  for  purposes  of  determining  the  characteristic  speeds,  which  are  the  eigenvalues  of  M~'  A.  After 
much  manipulation,  the  characteristic  polynomial  may  be  written  as  follows: 

2 


P{A)  =  [uf-  /i]2[  (u,  -  ][(«,  -  -*)-  -tt-t] 


Cl~fd f  +  1 


a\9f  + 1 


°2A  + 1 


<2  2  <9  +  1 J 

5  S 


with 


a2  = 


Ps  das 


and 


af  = 


Pfhf 

\t 

pfpfl 

hf\ 

\p  J  \ 

\r~Pf\ 

Mrh r\ 

-1) 

'P 

(39) 


(40) 


(41) 


Let  9r  -  — t - —  and  9  =  — r - - ,  where  VR  ,  and  V„  c  are  reference  velocities  for  the  fluid  and 

J  t/2  2  s  jr 2  2  7 


RJ 


Vh  a 


solids  phases.  With  these  definitions,  the  eigenvalues  may  be  found  as 

A  =  u/,u/  i[M/(l  +  M2,/)±^/(l-M2  /)2+4F/i2/],  \[u, (1  +  Ml ) ±  Ju, (1  - M2Rs )2  +  4v[s  ]  (42) 


where  the  reference  Mach  numbers  are  defined  as  M\  f  =  -^and  M\  s  =  — y- .  The  common  practice 

in  a  preconditioning  method  of  this  type  [21-23]  is  to  choose 
Vlf  =  min[max(w2  ,U2R),a2f] 

vls  =  min  [max(i/2  ,U\),a]] 

This  renders  the  characteristic  speeds  well  conditioned  at  all  flow  speeds.  Generalizing  this 
development  to  wave  propagation  in  the  direction  of  a  normal  vector  n  ,  we  have 


(43) 


/l  =  Vf  ■  n ,  Vf  •  n ,  -  [Vf  ■  n(l  +  M )  ±  ^Vf  •  n (1  -  M 2 4 )2  +  w[f  ], 


(44) 


K  ■  n, -IK -n(  1  +  M l ) ±  ^  •  "(1  - Ks  f  +  Ks  ] 


It  should  be  noted  that  one  can  consider  the  non-Newtonian  part  of  the  granular  stress  tensor  in  the 
eigenvalue  analysis.  In  this  case,  the  sound  speed  also  has  a  directional  dependence: 


dp > 


Ps  das 


1- 


sin(^)51>n,ny 


S  :  S  +  s / d2p  j 


(45) 


where  n,  are  the  elements  of  the  normal  vector  n  .  This  definition  may  be  used  in  the  finite-volume 
upwinding  methods  discussed  next  by  associating  the  normal  vector  with  that  pointing  outward  from  a 
cell  face. 


3.2.  Low  Diffusion  Flux-Splitting  Methods 

In  this  work,  we  extend  upwinding  techniques  previously  developed  for  isothermal  gas-solid  flows  in 
[7,8]  to  the  present  multi-phase  system.  Considering  again  one-dimensional  flow  for  simplicity,  we 
first  define  a  ‘numerical  speed  of  sound’  at  a  cell  interface: 


(46) 


V.l/2 


y«/(l  + 


1  +  M 


R.f 


'1/2 


where  the  'A  notation  indicates  the  evaluation  of  Eq.  (46)  using  cell-averaged  velocity  and  sound-speed 
data.  Other  quantities  needed  in  the  formulation  are  Mach  numbers  at  left  (L)  and  right  (R)  states: 

(47) 

af,V2 

and  Van  Leer  /  Liou  polynomials  in  Mach  number: 

M*  =Um±\M\), 


= 


[±i(M  ±  l)2,  \M\<\ 


(48) 


\M. 


(i)’ 


\M\>\ 


Given  these  basic  definitions,  the  convective  and  pressure  portions  of  the  interface  flux  are  expressed 
as: 

7”/,l/2  =  Pf.L&fJJ  P).Raf.R^  fin 


pp  -  P  q> 

rfM2  r/,l/2T 

with 

</>  =  [YJ,uf,H,]T, 

T  =  [0,0,l,0]r 

The  split  velocities  (7*  and  interface  pressure  PfU2  are  defined  as  follows: 

W  =  apvdKu  ~  M\n(\  1  )). 


(49) 


(50) 


2 , 


1/2 


U-  =  af  l/2[M~2)R  +  MU2(  1  +  ^  J^1  )], 


^■Pf.R^R.fM 


(51) 


^/,l/2  —  2  +  Pr)  +  2  P  Pr)  +  +  P  1) 

where  A p  =  pL  -  pR  ■  The  function  Mu2  is  given  as 

M,„  - -  M*u  -  M'„  +  M*„) 

and  the  Mach  number  functionals  used  in  the  interface  pressure  definition  are  given  as 


Pi  = 


U(l±M),  \M\<\ 


_l  ui 

M  '”(1)’ 


A/  |>  1 


(52) 


(53) 


Discretization  of  the  solids  continuity  term  and  the  interface  solids  pressure  follows  in  a  similar 
manner.  The  ‘numerical  speed  of  sound’  for  the  solids  phase  is  defined  as 


as,\/2  ~ 


JuliX-MlJ+AVl 


•  +<s 


(54) 
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and  the  interface  Mach  numbers  and  Mach-number  functionals  (Eqs.  (47),  (48),  and  (52))  arc 
redefined  using  the  solids  sound  speed  as  and  the  solids  velocity  us .  The  convective  flux  in  the  solids 
continuity  equation  and  the  interface  solids  pressure  then  become: 


^1/2  on,  =  Ps.Las.LU+  +  Ps.*as.RU~’ 

pP  -  D  IL» 

1  5,1/2  Vs, 1/2  T 

The  split  velocities  f/±and  interface  pressure  Ps  U2  are  re-defined  as: 


(55) 


(56) 


There  are  some  differences  in  the  split  velocity  and  interface  pressure  definitions  for  the  fluid  and 
solids  phases.  In  the  fluid-phase  forms,  the  pressure  diffusion  terms  in  the  continuity  equation  arc 
multiplied  by  af  so  that  in  the  limit  of  compaction  to  a  solid  state,  the  convective  fluid  flux  vanishes. 

The  solids  pressure  terms  naturally  vanish  in  the  limit  that  «vgoes  to  zero,  so  there  is  no  need  to 
multiply  the  pressure  diffusion  terms  in  the  solids  continuity  equation  again  by  as.  The  interface 
pressure  definition  for  the  fluid  phase  does  not  depend  on  the  volume  fraction,  as  the  entire  pressure 
gradient  will  vanish  as  a,  goes  to  zero.  For  the  solids  phase,  it  is  necessary  that  each  term  in  the 

interface  pressure  definition  vanish  as  a(goes  to  zero,  and  as  such,  the  velocity  diffusion  term  in  Eq. 
(56)  must  be  multiplied  by  as. 

3.3.  Time  Integration  Scheme 

The  three-phase  Navier-Stokes  system  described  above  is  implemented  into  a  version  of  REACTMB- 
MP,  a  general  purpose  Navier-Stokes  solver  for  multi-phase,  multi-component  flows  developed  at 
North  Carolina  State  University.  The  governing  equations  are  discretized  in  a  cell-centered  finite- 
volume  manner,  with  the  LDFSS  formulation  described  above  used  for  the  inviscid  fluxes.  Slope 
limited  total-variation  diminishing  (TVD)  variable-extrapolation  techniques  arc  used  to  extend  LDFSS 
to  second-order  accuracy  for  most  of  the  calculations.  One  set  of  calculations  for  water  entry  uses  the 
4th  order  accurate  Piecewise  Parabolic  Method  (PPM)  [24]  in  place  of  the  TVD  scheme.  Time 
integration  is  facilitated  by  a  planar  relaxation  sub-iteration  procedure  which  results  in  second-order 
temporal  accuracy  with  sufficient  sub-iteration  convergence.  For  steady-state  problems,  local  time- 
stepping  is  used  to  accelerate  convergence.  The  code  is  written  to  handle  simply-connected,  multi¬ 
block  meshes  and  is  parallelized  using  domain-decomposition  /  MPI  message-passing  methods.  To 
formulate  components  of  the  stress  tensor  at  cell  interfaces,  we  require  gradients  evaluated  at  the 
interfaces.  Given  a  cell  interface  i+1/2,  we  first  define  the  vector  pointing  from  cell-centers  i  to  i+l  as 


*  =  (*,+i  "  x,  )i  +  (yM  -  y,  )j  +  (zM  -  z, )k 

then  calculate  the  gradient  of  any  quantity  <j>  at  the  cell  interface  by 


(57) 


(58) 


r  r  r  r 
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where 


V^  =  |(V^+V4,)  (59) 

Gradients  at  cell  centers  are  calculated  using  standard  Green-Gauss  techniques.  The  splitting  of  the 
gradient  vector  into  a  component  aligned  with  f  and  one  perpendicular  to  f  helps  to  maintain  smooth 
solutions. 


3.4.  Arbitrary  Lagrangian-Eulerian  (ALE)  Formulation 


To  account  for  projectile  deceleration  due  to  hydrodynamic  forces,  we  have  adopted  an  Arbitrary 
Lagrangian-Eulerian  (ALE)  framework  [25-27].  In  this,  a  reference  frame  attached  to  the  projectile  is 
allowed  to  move  at  speeds  dictated  by  the  solution  of  Newton’s  laws  of  motion: 
dx„ 


—  —  u  „ 


dt 


m 


dup 

dt 


(60) 


-2* 


Integration  of  these  equations  at  each  time  step  provides  the  new  position  of  the  center  of  mass  of  the 
projectile  ( x  ),  referenced  to  a  datum  line  fixed  at  the  water-air  or  sand-air  interface,  and  the  new 


velocity  ( up  ),  which  is  used  as  a  fixed  grid  speed  in  the  ALE  formulation.  Velocity  components  arc 

measured  with  respect  to  this  reference  frame.  This  requires  some  modifications  to  the  interface  fluxes 
and  Jacobian  elements,  as  discussed  in  [26]  and  [27].  Though  the  ALE  formulation  generalizes  to 
three  dimensions,  only  rectilinear  deceleration  is  considered  in  this  paper. 


3.5.  Overset  Mesh  Framework 


Economical,  accurate  simulations  of  the  dynamic  entry  of  a  projectile  into  different  substrates  requires 
some  means  of  resolving  the  flowfield  around  the  body  while  accounting  for  the  long  distances 
(relative  to  the  penetrator  size)  that  must  be  traversed  before  the  projectile  stops.  Forces  and  moments 
generated  on  the  projectile  may  result  in  deviations  from  the  initial  path,  and  it  is  necessary  that  the 
computational  grid  follow  the  motion  of  the  projectile.  The  overset  mesh  approach  [28,9,10]  is  one  of 
the  more  powerful  techniques  available  to  handle  such  situations.  In  overset  meshing,  a  background 
mesh  with  coarser  spatial  resolution  is  constructed  for  the  entire  computational  domain.  Local  meshes 
which  may  move  relative  to  the  background  mesh  are  generated  around  the  projectile.  Flow  solutions 
within  the  background  and  local  meshes  are  time-advanced  in  tandem,  and  the  solutions  arc  transferred 
among  the  meshes  using  interpolation  procedures.  The  recent  development  of  automatic  C-language 
libraries  for  facilitating  this  information  transfer  [9,10]  has  greatly  simplified  the  addition  of  an  overset 
mesh  procedure  into  computational  codes.  The  SUGGAR  code  starts  with  an  assembly  of  overset 
meshes,  determines  which  meshes  lie  inside  one  another,  cuts  holes  in  the  overlying  meshes,  and 
determines  interpolation  stencils  for  transferring  information  from  one  mesh  to  another.  The  DiRTLIB 
routine  uses  information  from  SUGGAR  to  construct  local  communication  routines  that  facilitate  the 
information  transfer.  The  DiRTLIB  routines  link  with  the  main  CFD  code  in  a  straightforward 
fashion,  only  requiring  such  information  as  the  set  of  flow  variables  and  the  array  sizes. 
Communication  (through  MPI)  is  done  separately  from  the  normal  inter-block  communication  used  in 
the  code.  “Blanked  out”  regions  in  a  computational  domain,  corresponding  to  parts  of  the  flow  that  arc 
solved  on  other  meshes,  are  handled  using  an  IBLANK  array,  which  forces  the  code  not  to  update  the 
fluid  properties  when  IBLANK  is  one.  The  SUGGAR  code  is  written  to  be  used  with  meshes  that 
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generally  overlap  one  another,  as  opposed  to  the  simply  connected,  non-overlapping  meshes  used  in 
our  calculations.  A  significant  amount  of  time  was  spent  in  working  with  SUGGAR  to  enable  it  to 
correctly  respond  to  background  and  local  meshes  that  themselves  were  composed  of  non-overlapping 
meshes.  The  solution  involves  incorporating  details  of  the  block-to-block  connectivity  into 
SUGGAR’s  input  deck.  This  is  a  time-consuming  task  for  meshes  with  large  numbers  of  blocks,  and 
SUGGAR’s  format  does  not  account  for  complex  inter-connections  among  mesh  blocks  that  might 
occur  in  practice.  Future  work  will  require  the  development  of  software  that  automatically  performs  a 
translation  between  our  connectivity  format  and  that  required  by  SUGGAR.  Figure  3  presents  two 
overset  meshes,  one  of  which  uses  a  set  of  overlapping  meshes  for  the  background  mesh  and  the  other 
which  uses  a  set  of  simply-connected  meshes.  The  topology  of  the  composite  mesh  is  rendered 
reasonably  well  in  both  cases,  but  the  holes  are  cut  in  clearly  different  ways. 


3.6.  Sharp  Interface  Capturing 


To  account  for  wave  breakage  effects  during  water  penetration  and  for  precise  estimation  of  loads 
encountered  during  water  cavity  inception  and  collapse,  it  is  necessary  that  the  discretization  scheme 
capture  phase  interfaces  in  as  few  cells  as  possible.  This  can  be  accomplished  by  blending  the  baseline 
higher-order  scheme  (TVD  or  PPM)  with  a  bounded  downwind  differencing  scheme,  such  as  THINC 
[29]  or  CICSAM  [30],  THINC  (Tangent  of  Hyperbola  Interface  Capturing)  uses  the  tangent  hyperbola 
function  as  a  model  function  for  a  discontinuous  volume  (or  mass)  fraction  within  a  mesh  cell.  The 
average  volume  fraction  at  cell  i  is  defined  as 


a,  = 


Ajc, 


Hi  +  tanh(/7(  -l/2  -*,))]<&, 
J  ~  Ax. 


(61) 


where  £ ,  y ,  and  /?  arc  parameters  that  control  the  sharpness  of  the  approximation  and  its  connection 
with  cell-averaged  information  at  neighboring  cells.  The  location  of  the  center  of  the  distribution  x  is 

determined  by  enforcing  the  equality  in  Eq.  (61),  given  the  other  parameters.  Other  model  functions 
can  be  used  -  a  linear  distribution  may  be  represented  as 


a,  =y-  f  )-[\-y  +  2yx]dx. 
Ax.  J  ~ 


1  *,-i, 


(62) 


jc  =  max[0,min[l, 


x  -  x 


1-1/2  ~ 


Ax, 


-*,]] 


where  again,  x,  is  determined  by  enforcing  the  equality  in  Eq.  (62).  It  is  clear  that  this  framework  is 


directly  consistent  with  TVD  or  PPM  -type  interpolation  methods  in  that  left-  and  right-state 
information  may  be  obtained  easily  from  the  assumed  form  for  the  discontinuity  within  a  cell.  These 
techniques  tend  to  always  sharpen  gradients,  and  in  some  cases,  this  can  lead  to  unphysical  ‘stair- 
stepped’  interface  profiles  and  the  shedding  of  spurious  pockets  of  material  (known  as  ‘flotsam’  and 
‘jetsam’  in  the  literature).  To  avoid  this  behavior,  one  can  blend  the  baseline  TVD  or  PPM  scheme 
with  a  sharp-interface  scheme  using  a  linear  weighting  procedure.  In  the  CICSAM  scheme  [30],  the 
weighting  parameter  IT  at  a  cell  interface  is  a  function  of  the  angle  between  the  captured  phase 


Vcr  •  n 

interface  and  the  normal  vector  to  the  mesh-cell  interface  (IT  =  IT( - )),  and  if  the  angle 

’  I  Va 


approaches  ninety  degrees,  the  method  shifts  to  the  sharp-interface  capturing  strategy.  Our  tests  have 
shown  that  this  procedure  works  well  on  uniform  meshes,  but  on  the  highly  stretched  meshes  necessary 
to  capture  viscous  layer  development  on  the  projectile  surfaces,  the  method  leads  to  oscillations.  As 
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such,  wc  have  implemented  an  alternative  weighting  that  considers  not  only  the  direction  of  the  phase 
interface  but  also  the  magnitude  of  the  gradient: 

W  =  W{  ) ,  (63) 

\Va  n  |+- 
A 

where  s  /  A  is  a  cut-off  minimum  gradient.  An  advantage  of  this  form  is  that  the  sharp-interface 
contribution  will  be  diminished  as  the  gradient  becomes  small,  relative  to  the  cutoff,  even  if  the  phase 
interface  is  parallel  to  the  cell  interface.  This  should  allow  a  more  realistic  representation  of  two-phase 
mixing  processes  that  may  not  be  discontinuous  at  the  mesh  scale. 

3.7.  Grid  Arrangement 

Three  projectile  configurations  are  considered  in  this  work.  Closeup  views  of  some  of  the 
axisymmetric  grids  are  shown  in  Figure  4.  Configuration  I  (Figure  4a)  is  a  generic  dart  shape 
representative  of  those  used  in  beach  /  surf-zone  clearance  applications.  Configuration  II  (Figure  4b) 
corresponds  to  one  of  the  ballistic  shapes  (Round  10)  considered  in  an  early  investigation  of  projectile 
penetration  into  dry  sand.  [31,32]  Configuration  III  (Figure  4c)  was  used  in  sand-penetration  tests 
conducted  by  Lockheed  Martin.  [33]  Table  3  summarizes  the  total  number  of  interior  mesh  cells  and 
blocks  in  each  grid.  Also  listed  are  specifics  for  two  additional  renderings  of  Configuration  III  -  one 
three-dimensional  and  the  other  axisymmetric  with  overset  meshes  (see  Figure  3). 

Table  3:  Grid  Configurations  and  Data 


Configuration 

cells 

blocks 

I 

78973 

16 

II 

48948 

28 

III 

20947 

20 

III-C  (3-D) 

678034 

48 

III-D  (overset) 

34760  (projectile) 

25032  or  5376  (background) 

29 

4.  Results 

4.1.  Water  Entry 

The  first  set  of  results  focuses  strictly  on  the  entry  of  an  axisymmetric  projectile  into  water.  The 
granular-flow  model  is  not  utilized  in  these  calculations.  Figure  5  presents  drag  and  drag  coefficient 
versus  depth  for  dart  entry  (Configuration  I)  into  water.  A  dramatic  increase  in  drag  is  noted  as  the 
dart  first  enters  the  water  and  a  cavity  is  established.  From  this  point,  the  drag  decreases  steadily  as  a 
supercavity  is  formed.  Neaves  and  Edwards  [1]  simulated  this  case  using  a  supercavitation  code 
developed  at  NSWC-PC.  The  results  presented  are  in  good  agreement  with  the  earlier  ones.  The  use  of 
a  finer  grid  in  the  current  calculations  does  result  in  a  narrower  region  of  peak  drag,  as  the  water  /  air 
interface  is  captured  more  sharply. 

Figure  6  shows  snapshots  of  the  flowfield  about  the  nose  (top  sequence)  and  tail  (bottom  sequence)  for 
a  dart  (Configuration  I)  undergoing  deceleration  in  water.  The  right  side  of  each  figure  is  the  bulk  fluid 
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density  while  the  left  side  is  the  pressure.  Initially,  the  water  is  at  420  m/s  and  the  dart  is  completely 
contained  within  the  cavity.  As  the  dart  decelerates,  the  cavity  eventually  collapses  behind  the  dart  and 
moves  forward  along  the  body.  This  case  illustrates  the  ability  of  the  ALE  formulation  to  account  for 
this  time-dependent  dart  response.  Figures  7  and  8  present  drag  coefficient  versus  time  and  dart 
velocity  versus  depth.  Deviations  from  analytical  results  obtained  assuming  a  constant  drag  coefficient 
are  noted  as  the  dart  velocity  decelerates  below  80  m/s.  An  initial  decrease  in  drag  coefficient  occurs 
as  the  cavity  first  collapses  to  the  dart  trailing  edge  but  the  drag  coefficient  then  rises  as  the  back  of  the 
projectile  becomes  more  wetted. 

Figure  9  shows  an  example  of  the  use  of  sharp-interface  capturing  methods  (Section  3.6)  in  simulating 
vertical  entry  of  a  dart  into  water.  The  dart  geometry  and  initial  speed  are  that  of  Configuration  III 
and  the  refined  projectile  mesh  [Configuration  III-D  in  Table  3]  is  used  without  the  background  grid. 
The  ALE  formulation  is  not  used.  The  snapshots  compare,  from  left  to  right,  solutions  obtained  using 
the  PPM  scheme  without  interface-sharpening,  the  THINC  scheme  (Eq.  (61))  and  the  linear 
reconstruction  scheme  (Eq.  (62)).  The  top  set  of  snapshots  shows  the  liquid  water  mass  fraction,  with 
the  darker  contours  representing  the  water.  The  bottom  set  of  snapshots  shows  the  velocity  field 
induced  in  the  wake  of  the  projectile.  The  thin  lines  represent  the  10%,  50%,  and  90%  contours  of 
water  mass  fraction.  The  solutions  obtained  using  the  sharp-interface  capturing  schemes  show  more 
fine-scale  features,  including  secondary  instabilities  along  the  cavity  interface  that  result  from  vortex 
interactions.  The  appearance  of  these  features  may  be  due  to  the  fact  that  surface-tension  effects  arc 
much  better  resolved  when  the  sharp-interface  capturing  methods  are  used.  The  wake  of  the  projectile 
is  broader  and  more  dynamic  for  the  THINC  and  linear  reconstruction  methods. 

4.2.  Sand  Entry  -  Configuration  II 

Figures  10  and  11  present  results  for  steady  flow  of  a  sand  /air  mixture  over  the  Configuration  II 
projectile.  [32,33]  The  ALE  formulation  is  not  used  in  this  calculation,  solids-pressure  Model  I  is 
used,  and  the  calculation  is  conducted  in  a  non-time  accurate  manner  using  local  time  stepping. 
Conditions  correspond  to  those  encountered  just  after  entry  into  the  sand  bed  ( u  =  742  m/s,  as  = 

0.614).  A  steady  solution  is  sought  at  these  conditions  to  provide  an  initial  estimate  of  the  drag 
coefficient.  Figure  10  presents  solids  voidage  contours  for  two  values  of  the  assumed  loading  voidage 
aj0but  with  the  solids  shear  stress  tensor  omitted  ( <j>  in  Eq.  (15)  set  to  0).  The  sand  is  compacted  at  the 

nose  of  the  projectile  to  a  voidage  of  around  0.825  in  both  cases,  and  a  large  cavity  filled  with  air  is 
formed  as  the  granular  mixture  expands  around  the  blunt  nose  of  the  projectile.  A  shock  front  forms 
ahead  of  the  projectile  and  reflects  from  the  upper  wall  of  the  container.  This  front  is  associated  with 
the  onset  of  compaction,  which  raises  the  sound  speed  of  the  granular  material  to  levels  comparable  to 
the  fluid  velocity.  There  is  little  difference  in  the  solutions  for  the  two  choices  of art0 .  More  substantial 

differences  are  shown  in  Figure  1 1,  which  illustrates  the  effect  of  including  the  solids  shear  stresses 
(with  0=28  degrees).  The  bow  shock  angle  is  not  changed  substantially,  but  the  location  of  the 
air/sand  interface  shifts  toward  the  surface  when  shear  stresses  are  included.  Some  striations  in  the 
interface  downstream  of  the  end  of  the  projectile  are  observed.  These  tend  to  persist  over  the  iteration 
sequence,  precluding  good  steady-state  convergence.  The  inclusion  of  the  shear-stress  terms,  while  not 
destabilizing,  does  not  always  yield  smooth  solutions  for  volume  fraction  and  solids  pressure.  At  the 
nose,  the  solids  pressure  reaches  a  maximum  value  of  5.7e8  Pa  for  the  case  with  no  shear  stresses,  as 
compared  with  local  values  of  up  to  7.8e8  Pa  for  the  case  including  shear  stresses.  Figure  12  shows 
projectile  velocity  versus  depth  and  time  for  the  experiment.  This  data  can  be  regressed  to  yield  an 
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effective  drag  coefficient  Cd  of  2.02.  Table  4  indicates  that  the  calculations  predict  a  value  of  about 
half  of  the  experimental  value  and  that  the  inclusion  of  shear-stress  terms  slightly  increases  the  average 
drag  coefficient.  The  drag  coefficient  including  shear  stresses  does  not  stabilize  at  a  constant  value  but 
oscillates  between  values  of  about  1.13  and  1.35.  Due  to  the  formation  of  the  cavity,  the  drag  force  is 
due  almost  exclusively  to  the  solids  pressure  acting  on  the  projectile  nose,  implying  either  that  more 
compaction  should  occur  under  these  conditions  or  that  the  modeled  solids  pressure  is  too  low  for  the 
amount  of  compaction  that  is  predicted  to  occur. 

Table  4:  Drag  Predictions:  Configuration  II 


Configuration 

Measured  Cd 

Predicted  Cd 

II  ( aj0=0.614,  no  solids  shear  stress) 

2.02 

1.050 

II  ( aj0=0.64,  no  solids  shear  stress) 

2.02 

1.045 

II  (at0=0.614,  solids  shear  stress) 

2.02 

1.13-1.35 

4.3.  Sand  Entry  -  Configuration  III  (Initial  Results) 

Initial  results  obtained  for  the  Lockheed-Martin  dart  (Configuration  III)  arc  shown  in  Figures  13-15. 
Conditions  again  correspond  to  those  just  after  entry  into  the  sand  bed  (wp=  305  m/s,  as  =  0.56).  Solids 

voidage  contours  in  Figure  13  illustrate  the  effect  of  varying  the  voidage  value  (ars0)  at  which  the  sand 

is  assumed  to  begin  elastic/plastic  compression.  Solids  shear  stresses  are  omitted  in  this  case.  Regions 
of  strong  compression  of  the  sand,  followed  by  the  formation  of  a  cavity  filled  with  air,  arc  indicated 
for  both  cases  (av0  =  0.56  and  as0=  0.6).  The  onset  of  strong  compression  appears  to  be  dictated  by  the 

Mach  number  based  on  the  solids  sound  speed  (shown  in  Figure  14).  For  the  case  where  as0  =  0.56, 

corresponding  to  a  partially  loaded  initial  state,  the  solids  sound  speed  is  higher  in  the  granular  media 
and  increases  further  due  to  compaction.  The  locations  at  which  the  solids  velocity  becomes  lower 
than  the  solids  sound  speed  define  the  upstream  extent  of  the  bow  shock  located  in  front  of  the 
projectile.  The  solids  phase  is  initially  unloaded  when  as0  =  0.6,  the  sound  speed  is  extremely  small, 

and  some  compression  occurs  before  the  sonic  line  is  reached.  The  shock  front  is  located  closer  to  the 
body  and  curves  more  strongly  downstream  of  the  nose  region. 

Figure  15  shows  the  effect  of  including  the  solids  shear  stress  tensor  for  as0=  0.56.  The  dominant 

effect  is  to  force  the  sand/air  cavity  interface  closer  to  the  body.  This  leads  to  partial  “wetting"  of  the 
forebody  section  with  sand  and  thus  to  a  higher  drag  coefficient  (see  Table  5).  Figure  16  shows 
experimental  velocity  versus  depth  data  for  two  cases  involving  projectile  penetration  into  sand.  This 
data  has  been  used  to  extract  a  drag  coefficient  of  1 .03  for  the  initial  stages  of  entry.  Table  5  shows  that 
again,  the  numerical  approach  undcrpredicts  the  drag  coefficient.  If  sand  contacts  only  the  nose,  then 
the  predicted  drag  coefficient  is  around  0.21 .  The  partially  “wetted"  case  (including  shear  stress  terms) 
shown  in  Figure  15  yields  a  drag  coefficient  of  0.472,  but  the  same  case  with  as0=  0.6  (initially 
unloaded)  yields  a  drag  coefficient  of  only  0.264,  as  again,  a  full  cavity  is  formed. 


Table  5:  Drag  Predictions:  Configuration  III  (Initial  Results) 
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Configuration  III 

Measured  Cj 

Predicted  Cd 

ai0=0.56,  no  solids  shear  stress 

1.03 

0.220 

as0=0.6,  no  solids  shear  stress 

1.03 

0.213 

as0=  0.56,  solids  shear  stress 

1.03 

0.472 

as0  =0.6,  solids  shear  stress 

1.03 

0.264 

4.4.  Sand  Entry  -  Configuration  III  (Later  Results) 

Much  of  the  work  performed  after  the  publication  of  [34]  focused  on  trying  to  determine  the  causes  of 
the  underprediction  of  the  drag  coefficient  and  on  developing  methods  for  improving  the  results.  As  a 
first  step,  a  three-dimensional  mesh  was  generated  for  one  sixth  of  the  Lockheed-Martin  dart 
configuration  (Configuration  III-C  in  Table  3).  The  fins  were  included  in  the  computational  domain. 
Calculations  of  the  three-dimensional  flow  using  the  same  model  as  used  in  Sections  4.2  and  4.3  were 
performed  with  and  without  the  granular  shear-stress  model.  Figure  17  presents  a  snapshot  of  solids 
voidage  contours  for  these  cases.  The  fins  do  not  penetrate  the  sand  cavity  interface  in  either  case.  The 
rotational  flow  induced  by  the  fins  entrains  more  sand  into  the  wake  of  the  projectile,  and  some  sand  is 
collected  onto  the  flat  leading  edges  of  the  fins.  As  shown  in  Table  6,  the  inclusion  of  three- 
dimensionality  does  not  substantially  improve  predictions  of  the  drag  coefficient.  Because  of  the 
expense  of  the  calculations,  work  on  the  3-D  configuration  was  discontinued  in  order  to  focus  on  trade¬ 
off  studies  performed  on  the  axisymmetric  meshes. 

The  next  step  involved  re-fitting  the  solids  stress  model  based  on  Lockheed-Martin  sand-compaction 
data  (Model  II  in  Section  2.2),  which  was  not  available  prior  to  the  publication  of  [34],  This  change, 
along  with  the  addition  of  a  Newtonian  regularization  term  to  the  solids-stress  tensor  (Section  2.2)  and 
a  more  accurate  formulation  of  the  solids  stress  tensor  at  cell  interfaces,  yielded  improved  results  for 
the  drag  coefficient  (Table  6).  Figure  18  shows  contours  of  solids  volume  fraction  obtained  using  the 
improved  solids-stress  model.  The  inclusion  of  solids  shear  stress  effects  forces  the  sand  cavity 
interface  closer  to  the  body  surface.  This  leads  to  more  ‘wetting’  of  the  forcbody  section  with  sand  and 
a  higher  drag  coefficient.  Even  with  these  improvements,  it  became  clear  that  only  the  complete 
wetting  of  the  forebody  section  (up  until  the  maximum  cross-sectional  area)  would  lead  to  a  more 
correct  prediction  of  the  projectile  drag  coefficient,  at  least  under  axisymmetric,  steady  conditions. 

Table  6:  Drag  Predictions:  Configuration  III  (Later  Results) 


Configuration  III  (later  results) 

measured  Cd 

predicted  Cd 

as  o  =  0.56 ,  no  solids  shear  stress 

1.03 

0.220 

as  o  =  0.56 ,  no  solids  shear  stress  (3-D) 

1.03 

0.215 

aso  =  0.56 ,  with  solids  shear  stress  (3-D) 

1.03 

0.360 

as  o  =  0.56 ,  with  solids  shear  stress  (original  form 
presented  in  AIAA  Paper  2006-1289) 

1.03 

0.472 

as  o  =  0.56 ,  with  solids  shear  stress  (new  coding  of 

solids  stress  tensor  with  regularization  and  with  new  fit  to 
solids  pressure  (Figure  1)) 

1.03 

0.668  (average) 
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4.5.  Sand  Entry  -  Configuration  III-D 


The  latest  results  obtained  for  the  Lockheed-Martin  dart  configuration  have  used  a  refined  mesh 
(Configuration  III-D)  that  eliminates  the  hollowed-out  region  at  the  end  of  the  dart.  The  mesh  was 
generated  using  less-severe  stretching  ratios  and  more  mesh  points  per  block  in  the  hopes  that  the 
elements  of  the  granular  stress  tensor  would  be  resolved  more  correctly,  thus  providing  improved 
steady  state  convergence  and  more  accurate  results.  The  latest  runs  also  used  the  new  solids-pressure 
model  (Model  III  in  Section  2.2),  which  was  designed  to  produce  a  larger  increase  in  the  sound  speed 
as  one  approaches  a  loaded  state  while  maintaining  a  good  fit  to  the  Lockheed-Martin  sand-compaction 
data.  Several  parametric  studies  were  performed,  varying  such  factors  as  the  initial  solids  loading,  the 
inclusion  or  non-inclusion  of  the  direction-dependent  sound  speed  (Eq.  (45))  in  the  upwind  method,  the 
choice  of  the  constant  Cs  N  scaling  the  Newtonian  viscosity  component  of  Eq.  (15),  and  the  use  of  a 

‘thin  layer’  form  of  the  solids-stress  tensor  that  neglects  the  gradient  components  parallel  to  the  f 
direction.  Table  7  summarizes  some  of  the  results  of  these  parametric  studies. 

Table  7:  Parametric  Studies  Performed  on  Configuration  III-D 


Case  number 
(Configuration  III-D) 

«*,<> 

CS,N 

directional 
sound  speed 

thin  layer  form 

predicted  Cd  (range) 

1 

0.56 

0.00 

no 

no 

0.529-0.543 

2 

0.56 

0.00 

yes 

no 

0.465-0.469 

3 

0.56 

0.01 

no 

no 

0.557-0.531 

4 

0.56 

0.10 

no 

no 

0.578-0.586 

5 

0.56 

1.00 

no 

no 

0.747-0.853 

6 

0.58 

1.00 

no 

yes 

0.740-0.810 

Many  other  combinations  of  parameters  were  tried.  In  general,  increasing  the  initial  value  of  as  toward 
the  fitted  value  for  the  onset  of  compaction  ( ocs,=  0.580769002438)  led  to  increases  in  the  drag 
coefficient,  as  did  increases  in  Cs  N  .  In  some  cases,  however,  this  response  was  accompanied  by  the 

destabilization  of  the  cavity  wake,  initiated  by  large  oscillations  in  the  vicinity  of  the  cavity  inception 
point.  Even  when  destabilization  did  not  occur,  convergence  was  very  poor.  A  significant  effort  was 
undertaken  to  understand  the  causes  of  convergence  degradation.  It  was  found  that  it  was  a  result, 
primarily,  of  oscillations  in  the  solids  properties  in  the  vicinity  of  the  projectile  nose  but  that  large  drag 
terms,  caused  by  a  wide  disparity  in  solids  and  fluid  velocities  in  the  projectile  wake,  also  contributed 
to  the  response.  Figure  19  shows  a  close-up  view  of  solids  volume  fraction  and  pressure  in  the 
vicinity  of  the  projectile  nose  for  Case  5.  The  formation  of  the  cavity  is  delayed  until  the  maximum- 
area  point,  and  the  entire  forebody  of  the  projectile  is  ‘wetted’  by  the  sand. 

4.6.  Overset  Mesh  Results  -  Configuration  III-D 

As  mentioned  in  Section  3.5,  we  were  able  to  embed  the  computational  model  into  an  overset-mesh 
framework  using  the  SUGGAR/DiRTLIB  [9,10]  utilities.  We  exercised  this  capability  for  a  few  cases, 
focusing  primarily  on  ensuring  that  the  procedures  provided  the  correct  information  transfer  between 
the  projectile  meshes  and  the  background  meshes.  Figure  20  shows  results  obtained  using  SUGGAR  / 
DiRTLIB  for  a  case  where  the  background  meshes  were  allowed  to  overlap  one  another  (bottom)  and 
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for  a  case  when  the  background  meshes  do  not  overlap  (top).  The  solutions  correspond  to  the  meshes 
shown  in  Figure  3.  The  variable  plotted  is  the  transverse  component  of  the  solids  velocity.  The  cases 
were  run  with  different  parameters,  and  close  agreement  between  the  solutions  is  not  to  be  expected. 
The  contours  indicate  that  information  is  being  passed  correctly  from  one  block  set  to  another.  The 
resolution  differences  between  the  background  and  projectile  meshes  lead  to  some  smearing  of  the 
waves. 

5.  Conclusions 

The  development  of  a  comprehensive  three-phase  model  for  simulating  projectile  penetration  into 
water  and  sand  has  been  outlined  herein.  The  approach  combines  low-diffusion  upwinding  techniques 
valid  for  multi-phase  flows  with  time-derivative  preconditioning  methods  to  ensure  efficient,  accurate 
time  evolution  at  all  speeds.  The  computational  model  has  been  embedded  into  an  overset-mesh 
framework,  which  renders  it  suitable  for  eventual  simulations  of  six-degree-of-freedom  maneuver. 
Applications  to  supercavitating  projectile  entry  into  water  and  subsequent  deceleration  show  good 
agreement  with  prior  work  and  with  experimental  trends.  Applications  to  situations  representative  of 
initial  projectile  entry  into  dry  sand  indicate  that  the  current  approach  underpredicts  measured  drag 
coefficients  by  more  than  75%  in  some  cases.  The  addition  of  a  solids  pressure  model  that  provides  a 
sharper  transition  to  a  loaded  state  upon  initial  compaction  and  the  inclusion  of  an  ad-hoc  Newtonian 
contribution  to  the  baseline  Mohr-Coulomb  granular  stress  tensor  results  in  better  agreement  with 
experimental  data,  with  the  best  predictions  for  the  drag  coefficient  being  -20%  lower  than  the 
measurements.  These  results  indicate  that  the  current  models  for  solids-stress  (frictional)  effects  need 
to  be  improved.  The  assumption  of  axisymmetric  flow  may  also  be  overly  restrictive,  as  it  is  known 
that  the  motion  such  projectiles  may  exhibit  sinuous  and  rotational  modes  accompanied  by  ‘tail  slap’ 
on  the  walls  of  the  cavity.  Complete  six-degree  of  freedom  simulations  of  dynamic  entry  would 
therefore  be  required  for  a  more  complete  analysis. 
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solids  pressure  (Pa) 


a.)  Model  1  b.)  Model  II 

Figure  1 :  Solids  pressure  versus  volume  fraction  for  Cooper-Eaton  Models  I  and  II 


Figure  3:  Solids  pressure  versus  volume  fraction  for  re-fitted 
Model  III 
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(UI)A 


Figure  3:  Composite  grids  generated  by  SUGGAR  for  overlapping  and 
non-overlapping  background  meshes 


Figure  4:  Grid  details  for  projectile  configurations  (grid  reflected  about 
centerline  for  clarity) 
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Figure  5:  Drag  and  drag  coefficient  versus 
time  for  dart  entry  into  water 
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Figure  6:  Snapshot  of  pressure  and  density  contours:  dart  entry  into  water 
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Figure  7:  Drag  coefficient  versus  time:  dart 
deceleration  in  water 
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Figure  8:  Velocity  versus  depth:  dart 
deceleration  in  water 


Figure  9:  Predictions  of  vertical  water  entry  using  various  interface-sharpening 
schemes  (from  left  to  right:  PPM,  PPM  with  THINC,  PPM  with  linear  reconstruction; 
top:  water  mass  fraction;  bottom:  vertical  velocity) 
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Figure  10:  Solids  voidage  contours  -  effect  of  asQ :  Configuration  II 
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Figure  1 1 :  Solids  voidage  contours  -  effect  of  solids  shear  stresses: 
Configuration  II 
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(Ui)  A 


■  Exper  Penetration  Data  (round  10) 

- Anal  (Cd  =  2.02) Time  vs  Depth 

-  Anal  (Cd  *  2.02)  Vel  vs  Depth 


Depth  (m) 

Figure  12:  Velocity  versus  depth  and  time  (experiment): 
Configuration  II 


Figure  13:  Solids  voidage  contours  -  effect  of  as0 :  Configuration  III 
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Figure  14:  Solids  sound  speed  contours:  Configuration  III 
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Figure  15:  Solids  voidage  contours  -  effect  of  solids  shear  stresses: 
Configuration  III 
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Figure  16:  Velocity  versus  depth  (experiment):  Configuration  III 


Figure  17:  Solids  volume  fraction  contours  for  3-D  simulation  of 
dry  sand  flow  over  the  Lockheed-Martin  dart  (including  fins)  Left 
-  no  solids  shear  stress;  right  -  with  solids  shear  stress 
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Figure  19:  Solids  volume  fraction  and  solids  pressure  contours  using  improved 
solids-stress  models. 
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Non-overlapping  background 
mesh 


Figure  20:  Solids  vertical  velocity  predictions  on  two  different 
composite  meshes 
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